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We present radiation transfer models that demonstrate that reflected light levels 
from three dimensional (3D) exoplanetary atmospheres can be more than 50% lower 
than those predicted by models of homogeneous or smooth atmospheres. Compared to 
smooth models, 3D atmospheres enable starlight to penetrate to larger depths result- 
ing in a decreased probability for the photons to scatter back out of the atmosphere 
before being absorbed. The increased depth of penetration of starlight in a 3D medium 
is a well known result from theoretical studies of molecular clouds and planetary at- 
mospheres. For the first time we study the reflectivity of 3D atmospheres as a possible 
explanation for the apparent low geometric albedos inferred for extrasolar planetary 
atmospheres. Our models indicate that 3D atmospheric structure may be an impor- 
tant contributing factor to the non-detections of scattered light from exoplanetary 
atmospheres. We investigate the self-shadowing radiation transfer effects of patchy 
cloud cover in 3D scattered light simulations of the atmosphere of HD209458b. We 
find that, for a generic planet, geometric albedos can be as high as 0.45 in some limited 
situations, but that in general the geometric albedo is much lower. We conclude with 
some explanations on why extrasolar planets are likely dark at optical wavelengths. 

Key words: radiative transfer, planetary systems, planets and satellites: individual: 
HD 209458b, atmospheric effects 



1 INTRODUCTION 

Since the discovery of the planet orbiting 51 Peg over a 
decade ago ( [Mayor &: Queloz|1995[ ), the detection and char- 
acterization of extrasolar planets has become an area of in- 
tense study. To date, over two hundred and seventy extraso- 
lar planets have been discovered primarily by the measure- 
ment of the doppler shifting of stellar spectral lines due to 
the motion of the star around the centre of mass of the exo- 
planetary system. Other indirect techniques that have suc- 
cesfuUy detected extrasolar planets include transit searches 
and microlensing (e.g. Pont (20061; Beaulieu et al. ('2006)). 



fleeted starlight from exoplanetary atmospheres is therefore 



a realistic goal (e.g., Seager et al. (2000l). However, despite 
intense effort observers have been unable to detect stellar 
light reflected from extrasolar planetary atmospheres. |Col-| 



lier Cameron et al. ( 1999 1 and Charbonneau et al. ( 1999h de- 



veloped techniques that exploit the doppler effect to separate 
the spectral lines of an extrasolar planet from the lines of the 
parent star, but both have been unable to definitively detect 
the faint signature. A detection of reflected light from the 



planet orbiting Tau Bootis by Collier Cameron et al 



was unverifiable in the following observing season 



The first detection of direct light from an extrasolar planet 
was the infrared emission due to reprocessed starlight in the 
warm atmospheres of the planets orbiting TrES-1 and HD 
209458 ( [Charbonneau et al. 12005) [Deminget al.|2005[ ). Addi- 
tionally, a handful of planets may have been detected using 
direct imaging ( Neuhaeuser|2005"| 

Many extrasolar planets orbit within a few astronom- 
ical units of their parent star. The direct detection of re- 



(19991 



Collier 



Cameron et al. 20041. Subsequent work has set an upper 



limit of 0.39 for the planet's geometric albedo, (Ag), as- 



suming a grey albedo and Venus-like phase function (Leigh 
et al.|2003K Using similar methods the group found the ge- 



ometric albedo of the companion to HD 75289 has an upper 
limit of Ag < 0.16 ( Leigh et al.|2003 |. Analysis of data from 
the MOST satellite has set an optical (400-70 nm) limit of 
Ag < 0.17 on HD 209458b ( [Rowe et al.p007 l. 



B. Hood et al. 



These non-detections of reflected light and upper Hmits 
on the geometric albedo imply that compared to solar sys- 
tem planets, the extrasolar planets are very dark at optical 
and near-infrared wavelengths. In comparison, the geometric 
albedos of solar system gas giant planets range from 0.41 to 
0.52 ( Cox||2000 \ . Particles in the exoplanetary atmospheres 
could have very low scattering albedos or the particles could 
have scattering phase functions that are very forward throw- 
ing so that starlight is directed deep into the atmosphere 
from where it has a small probability of scattering back out. 



Sudarsky et al. ( 2000 \ explored parameters which can influ- 



ence the geometric albedo of an extrasolar planet, with a 
one-dimensional radiation transfer model. The focus of our 
paper is the effect of three-dimensional atmospheric struc- 
ture on the radiation transfer of starlight and hence the over- 
all reflectivity of exoplanetary atmospheres. 

Previous studies have shown the importance of three- 
dimensional radiation transfer effects on the penetration of 
starlight into dark clouds ( Boisse||1990 l, the analysis of re- 
flection nebulae ( Mathis et al.|2002 l, and the infrared spec- 
trum from dusty ultracompact H ll regions (Indebetouw 
et al.|[2006 |. In the planetary atmosphere community much 
effort has been devoted to developing three-dimensional ra- 
diation transfer techniques to study the penetration of Solar 



radiation through clouds in the Earth's atmosphere ( Caha- 
lan et al.|2005 |. Studies of clouds in the Earth's atmosphere 
have shown a fractal structure Lovejoy ( 1982 1. Additionally, 



Kuchinke et al. ( 2004 1 showed that a three-dimesional Monte 



Carlo radiation transfer simulation in fractal clouds results 
in a more accurate representation of radiation penetration 
in real clouds than a plane-parallel homogenous simulation. 
Though we lack the large amount of data of atmospheric 
scientists studying the Earth, this research is very relevant 
to extrasolar planetary atmospheric simulations. 

The work we present in this paper is motivated by 
the non-detections of reflected light from exoplanets and 
the clear three-dimensional nature of clouds in the atmo- 
spheres of the Earth and other planets in the Solar Sys- 
tem. We present three-dimensionalj radiation transfer sim- 
ulations of the reflected light from inhomogeneous exoplan- 
etary atmospheres. We have not computed the 3D temper- 
ature and pressure structure of an atmosphere, rather to 
investigate 3D radiation transfer effects on reflected light 
we take ID atmospheres and convert them to 3D via a hi- 
erarchical clumping algorithm. This approach allows us to 
identify and characterize the different reflectivities of uni- 
form and 3D atmospheres. We then extend our analysis to 
study a more detailed atmospheric model incorporating mul- 
tiple species of scatterers, cloud decks, and vertical density 
gradients within the atmosphere. Although we do not solve 
for the 3D temperature, pressure, and density structure, we 
attempt to build a more realistic simulation and will focus 
specifically on current atmospheric models for the extraso- 
lar planet HD 209458b. We have simulated the atmosphere 
and found that in all cases the geometric albedos for the 
simulations of extrasolar planet HD 209458b are lower, and 
in many cases much lower, than that of Jupiter, suggesting 



^ We use 3D to mean that we can vary the density and charac- 
teristics of our model for each point in our x,y, and z dimensional 
grid, and introduce structure that varies in every direction. 



radiation transfer effects in a 3D atmosphere might help ex- 
plain why to date, reflected light has not been detected from 
an extrasolar planet. 

In !j2]we describe our numerical radiation transfer code, 
the algorithm we use for generating a hierarchically clumped 
3D atmosphere, and the adopted scattering properties of 
particles within the atmosphere. Results of our radiation 
transfer simulations for a wide range of atmospheric opti- 
cal depth, porosity, and particle scattering properties are 
presented in !|3]and §3. 4| summarizes our generalized results 
with respect to the scattered light levels. !|4] presents 3D 
radiation transfer models for the scattered light from the 
atmosphere of HD 209458b and in 35] we summarize our 
findings. 



2 MODEL INGREDIENTS 

This section describes the radiation transfer technique, at- 
mospheric geometry, and particle scattering properties we 
adopt for our 3D radiation transfer simulations. 



2.1 Monte Carlo Radiation Transfer Code 

For our refiected light simulations we use a 3D Monte Carlo 



scattering code, described in Wood & Reynolds ( 1999 1. This 



code simulates radiation transfer through a 3D linear Carte- 
sian grid onto which the density structure is discretized. 
The code accurately treats multiple, anisotropic scattering 
for any analytic or tabulated angular scattering phase func- 
tion. We do not consider the reprocessing of absorbed stellar 
photons. This assumes that all photons that are absorbed by 
the planetary atmosphere are re-emitted at wavelengths that 
are well separated from the incident wavelength. Therefore 
our simulations will be accurate at optical and near-infrared 
wavelengths where refiected light dominates the planetary 
spectrum and thermal emission from the planetary atmo- 
sphere is negligible. To construct a simulation for a refiected 
spectrum using our technique we would run our code multi- 
ple times (wavelength by wavelength) and change the scat- 
tering properties of the medium to be those appropriate for 
each wavelength in the incident spectrum. 

The density grid we employ for the simulations in this 
paper has 200 cubical cells on a side, for a total of 8 x 10^ 
cells. Our atmosphere is semi-infinite, so that photon pack- 
ets exiting the grid through x or y faces, re-enter the grid on 
the opposite face. Thus photons can only escape through the 
upper and lower z boundaries of the grid. See § |2.4| below 
for our treatment of different lower z boundary conditions. 
Because of our semi-infinite treatment of the atmospheric 
simulation, we are unable to reproduce the small amount of 
light which might otherwise pass through a tenuous atmo- 
sphere when the planet is directly between the star and the 
observer (at a 180°phase angle). Since we are considering 
reflected light and not transmission spectra (see Seager & 
Sasselov] pOOO ) ) , this is not an obstacle to our approach. 

The grid is illuminated from the outside with photons 
entering the upper z face at specified angles representing the 
incident angle of stellar photons on the atmosphere. When 
photons exit the grid through either the upper or lower z 
faces, they are binned in solid angle according to their di- 
rection of travel. The output of a simulation is the direction 
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Figure 1. Reflected light from one of our semi-inflnite simula- 
tions. The density grid is uniformly illuminated by plane parallel 
rays at an angle of 20° to the surface. The reflected light is plotted 
as a function of longitude (rp) and latitude {fi = cos 6) around the 
grid, and is normalized to isotropic, non-absorptive scattering. 



dependent reflected light levels in 400 evenly spaced solid 
angle bins (20 divisions in cos 6 and 20 divisions in (f)). We 
ran simulations for each realization of the atmosphere geom- 
etry with incident angles in the range 0° to 85° (measured 
from the normal to the upper z face of the simulation grid) 
in 5° increments. The results of those simulations are the 
flux emitted in each of the 400 solid angle bins, and can be 
presented graphically, as in flgure [T] 

Following other theoretical simulations of reflected light 
(e.g., Sudarsky et al. (2005')), the results of our reflected light 
simulations from the semi-inflnite models are tiled to gen- 
erate the total reflected light from the planet. In the semi- 
inflnite simulation we specify the incident angle of radiation 
and the output is the angle dependent reflectivity. Using 
the results from the semi-inflnite simulations we may then 
specify an incident angle of radiation, an emergent angle of 
scattered radiation and flnd the reflected light level com- 
pared to a Lamb ertiarp] scattering surface. Using equation 
9.8 fromlSobolevI (pTsf 



Hia) 



r/2 



cos{a—uj) cosu!<kj 



it/2 



r/2 



p(77,C,0)cos ipdip{l) 



we determine H{a), the radiation flux emerging from a 
planet with respect to the planetary phase a. In this equa- 
tion 7] is the cosine of the angle of emergence with respect 
to the outward normal, ^ is the cosine of the angle of in- 
cidence with respect to the outward normal, is the az- 



^ A Lambertian surface is completely non-absorbing and scatters 
an equal intensity in all directions. 



imuthal angle between those two angles. The planetocentric 
coordinates of latitude and longitude are represented by ip 
and Lu. Finally, p is the output from our Monte Carlo sim- 
ulations and it is the reflection coefflcient, or the flux scat- 
tered in a particular direction from the surface compared to 
an isotropically scattering surface. Our treatment of tiling 
results from many semi-inflnite simulations follows that of 



Sudarsky et al. (20051, but we use Monte Carlo techniques 



for the 3D radiation transfer. 



2.2 3D Atmospheric Density Structure 

The goal of this paper is to investigate the effects of 3D 
density structures and 3D radiation transfer on the reflec- 
tivity of extrasolar planets. Therefore we begin with the sim- 
plest geometry of a uniform density atmosphere (no vertical 
stratiflcation) and compare the resulting reflectivity with 3D 
atmospheres of the same total mass. As mentioned in the 
introduction, we do not compute the 3D density, pressure, 
and temperature structure of an atmosphere, but convert 
smooth ID atmospheric densities to 3D using the hierarchi- 



cal clumping algorithm of Elmcgreeii|(1997l as described in 



Mathis et al. ( 2002 \ and modified in I Wood et al. ( 2005 I 



ien| 
ml 



The hierarchical clumping algorithm randomly casts A'^i 
points in 3D space at the first level. At the second level 
A*' points are cast around each of the TVi points from the 
first level, but with a separation in the range ±A'^~^'/2, 
where H is the hierarchical level and / = log N/ log A is 
the fractal dimension of the hierarchical structure. At all 
subsequent levels a further A'' points are cast around each of 
the points cast at the previous level, again with separation in 
the range ±A'^~^Y2- The density in our Cartesian grid is 
proportional to the number of points in each cell that were 
cast at the final hierarchical level. See Figure 5 in |Wood| 



et al. ( 2005 1 for a graphical representation of the algorithm 



and how the density in each cell is determined. 

In the simulations presented here we use five hierarchi- 
cal levels, TV — 32, and fractal dimension / — 2.6. This 
value of the fractal dimension approximately corresponds 
to a projected two-dimensional area-perimeter dimension 
of 1.36, appropriate for clouds in the Earth's atmosphere 
( Lovejoy|1982| . 

We convert the uniform density to 3D using this algo- 
rithm and leave a fraction of the mass, /smooth, in a smooth 
density component, the remainder being distributed accord- 
ing to the hierarchical algorithm. The minimum vertical op- 
tical depth for a 3D model is r^^ = /smoothTz , where t^ is 
the optical depth of the smooth model and /smooth is the 
mass fraction smoothly distributed in the 3D model. This 
optical depth occurs for sightlines along which there is no 
material distributed by the hierarchical clumping algorithm. 

To investigate different porosity levels of the atmo- 



sphere we adopt the approach of Wood et al. ( 2005 1 and 



vary the number of points, A'^i, cast at the first hierarchical 
level. For small values of A'^i the porosity will be high and as 
Ni increases the density approaches a more uniform distri- 
bution. Specifically we investigate atmosphere models with 
A'^i = 2, 8, 32, and 128. These values correspond to poros- 
ity values of around 99%, 95%, 80%, and 45% respectively, 
where we define porosity as the fraction of the volume occu- 
pied by the smooth density component. We also determine 
the covering factor of the hierarchically distributed clumps 
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in the grid by creating column density maps viewed along 
the 2 axis. For the TVi values used in the paper the cover- 
ing fractions are approximately 15%, 45%, 90%, and 100% 
respectively. These values of the porosity and cloud cover- 
age do not change significantly (within a few percent) for 
different random castings of the first Ni points in the hier- 
archical clumping algorithm. For the rest of the paper we 
refer to different 3D simulations according to their porosity 
or covering factors. In the paper we refer our 3D results to 
"the equivalent smooth model," which means the smooth 
model with a particular optical depth (i.e., uniform density 
or mass) which is then converted to 3D using the hierarchical 
clumping algorithm. 

Figures [2] and [3] show examples of the hierarchical den- 
sity grid. Figure [2rshows column density maps for the four 
different porosity values we investigated. The lowest column 
density in each map arises when there are no hierarchically 
distributed clumps along the line of sight, and so in these 
maps black is not zero density, but the column density due 
to the smooth component only in the density grid. Figure [3] 
shows a voxel projection of the density grid for a poros- 
ity value of 80%. Note that our implementation of a semi- 
infinite grid for the radiation transfer simulation is equiva- 
lent to an infinite sequence in x and y of such density cubes. 

We have the ability to change our grid size in order to 
investigate finer or coarser resolution, but with current set- 
tings our individual simulations represent approximately a 
1200 X 1200 X 1200 km cube on the surface of a Jupiter-sized 
planet. When we simulate HD209458 in § [4l we adjust our 
simulation to investigate a deeper portion of the atmosphere 
with a grid of 8000 x 1600 x 1600 km. 



2.3 Particle Scattering Properties 

The Monte Carlo code simulates the scattering of stellar 
photons off particles in the planetary atmosphere. The code 
can treat analytic or tabulated scattering phase functions. 
In our initial investigation we represent the angular shape 
of scattering with the Henyey-Greenstein phase function 



Henyey & Greenstein ( 1941 1 characterized by the single pa- 
rameter, g, 



HG{v) = 



1 



2(l + p2-2gcosu)3/2 



(2) 



where v is the scattering angle. The scattering asymmetry 
parameter controls the shape of the phase function, with g — 
giving isotropic scattering, positive g for forward throwing, 
and negative g for a phase function dominated by back- 
scattering. 

The scattering albedo a, determines the fate of photons 
at their interaction locations in the atmosphere. Absorbed 
photons are terminated and assumed to be reprocessed to 
long wavelengths and do not contribute any flux at the sim- 
ulation wavelength. 

Typical values for the albedo and scattering asymmetry 
parameters in the optical are a ~ 0.5 and g ~ 0.5, based 
on Mie theory calculations for grain types assumed to be 
present in extrasolar planetary atmospheres. Enstatite, iron, 
and corundum grains are expected in close-in giant planet 
atmospheres (Seager et al.|2000l. 




Figure 2. Vertical column density maps for four different poros- 
ity levels produced by our hierarchical clumping algorithm. We 
define porosity to be the fraction of our simulation grid that con- 
tains only the smooth density component. The figures show poros- 
ity levels of 99% (top left), 95% (top right), 80% (bottom left), 
and 45% (bottom right). 



2.4 Lower Boundary Conditions 

Simulated photons are injected into the density grid as 
described above and may exit the grid through the up- 
per z face in which case they are placed into direction-of- 
observation solid angle bins. For photons that exit the simu- 
lation through the lower z face we investigate two scenarios, 
absorption or Lambert reflection. In the case of absorption, 
this may represent the absorption of photons by an optically 
thick layer below the cloud deck. The second case is equiva- 
lent to photons being isotropically scattered off a planetary 
surface; the simulation returns the light curve expected for a 
Lambert sphere when the optical depth of the atmosphere is 
zero (see 



2.6 below). These two boundary conditions have 



the largest influence on the reflected light for atmospheres 
where photons penetrate to the lower z boundary, such as 
optically thin atmospheres and those with very low porosity. 



2.5 Albedo definitions 

We use two deflnitions of albedo throughout this paper. 
The first albedo, called single-scattering albedo, refers to 
the albedo of the particles in the atmosphere. We assign 
an albedo to each scattering species. For instance, if we as- 
signed an albedo of 0.6 to a scattering species representing 
iron, then approximately 60% of the Monte Carlo photons 
interacting with iron will be scattered, and 40% absorbed. 
The second albedo we use is the geometric albedo, Ag, 
which describes the albedo of the entire planet. The geomet- 
ric albedo is defined relative to the amount of light refiected 
from a flat, Lambertian disk with the same radius as the 
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Figure 3. A voxel projection of the 3D density distribution for 
a density grid with porosity of 80%, corresponding to the lower 
left panel in Fig. [2] 



planet, at the same distance from the star as the planet. 
Comparing geometric albedo to the Lambert result is only 
strictly correct when the planet is at a phase angle of zero 
degrees. However, the reflected hght at the other phase an- 
gles can be compared to this zero phase geometric albedo 
and we present this in our figures below. 



2.6 Code Testing 

We have tested our Monte Carlo scattering codes against an- 
alytic solutions and other independently developed radiation 
transfer codes. For scattering simulations, Chandrasekhar's 
solution for a plane parallel Rayleigh scattering atmosphere 
illuminated from below is one of the standard benchmarks. 
Our codes accurately reproduce Chandrasekhar's results for 
the angular dependence of the emergent radiation (intensity 
and polarization) from the top of the atmosphere (e.g., see 
Wood et al. ( 1996 1). In this paper we wish to tile results from 
a semi-infinite Monte Carlo simulation to generate a phase 
function for the reflected hght from the planet as a whole. 
To test that we correctly implement the tiling described in § 
|2.1| we set up a simulation to represent a Lambert sphere. In 
this simulation the atmosphere has zero optical depth and 
photons at the lower boundary reflect isotropically. 

Our tiling procedure reproduces the Lambert results to 
within 1% at small phase angles, and to much better accu- 
racy at larger phases. These small discrepancies are due to 
the coarseness of our numerical integration. Increasing the 
angular resolution of the integration by a factor of ten repro- 
duces the Lambert phase function to within 0.01%. However, 
the small increase in accuracy does not justify the significant 
increase in computational time, so the results presented in 
the rest of the paper use the lower angular resolution. 



2.7 Varying the cloud distribution for the same 
porosity 

To set up the 3D density structures we use an initial random 
number as a seed for the hierarchical clumping algorithm. 
Different random seeds produce different arrangements of 
clouds for the same porosity value and most likely will yield 
different levels of reflected light. We have investigated the 
variations in reflectivity arising from different 3D structures 
with the same porosity levels and find the largest differences 
(up to 20%) occur for simulations with very high porosity. 
In these situations there are few cloud complexes in the sim- 
ulation grid and if clouds lie above one another the cover- 
ing factors can be very different for the same porosity level. 
Smoother atmospheres exhibit much smaller variations (typ- 
ically 5%) among different cloud geometries with the same 
porosity. 



3 3D SCATTERED LIGHT SIMULATIONS 
3.1 Fiducial Smooth Model 

We have chosen a set of uniform density fiducial models 
against which we compare reflected light levels with those 
from the 3D models. The flducial models have a uniform 
density, phase function asymmetry parameter g — 0.5, par- 
ticle scattering albedo a = 0.5, and total vertical optical 
depths through the grid of r^ = 0.1, 1, 10, and 100. The op- 
tical depths and scattering parameters of the fiducial models 
are similar to those expected in extrasolar planetary atmo- 



spheres ( Seager et al. 2000 1 and provide a benchmark against 
which we can explore the effects of 3D density structures. 
The lower boundary condition for these fiducial models are 
that all photons at the lower boundary are absorbed. 

Figure |4] shows the derived geometric albedos as a func- 
tion of phase angle for the tiled simulations for these four 
smooth models. An immediate result from the fiducial mod- 
els is that the reflectivity and hence geometric albedo does 
not change for optical depths greater than r^ = 10. At high 
optical depths the reflectance is due to scattered photons 
that only penetrate to a small physical depth in the atmo- 
sphere. The geometric albedo is dominated by photons that 
are scattered off the top of the optically thick atmosphere. 



3.2 Parameter Exploration 

The main focus of our paper is to investigate the effects on 
reflected light levels of 3D radiation transfer in clumpy atmo- 
spheres. However, as has been reported in other theoretical 
reflected light studies of ID atmospheres, many parameters 



can affect the overall reflected light levels (Sudarsky et al. 
|2000| [Barman et aL||2005| [Sudarsky et al.||2005| ). As~de^ 
scribed in detail in Hood (20071, we have explored a very 



wide parameter space and investigated the effects of varying 
particle albedo, phase function, and lower boundary condi- 
tions. In this section, we first briefiy describe these results 
before presenting figures for our main study: the effects of 
atmospheric porosity on refiected fight levels from extrasolar 
planets. 
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Figure 4. Effect of r^ on reflectivity for our set of fiducial models 
at a range of optical depths. All models have single particle albedo 
a = 0.5, an absorbing lower boundary, and an atmosphere with a 
uniform density structure. For r > 10 the results do not change. 



3.2.1 Particle albedo and phase function 

One of the most significant parameters set within the sim- 
ulation is the albedo of the scattering particles in the at- 
mosphere. This was set to a = 0.5 in the fiducial models 
in Figure [4] and produced a geometric albedo for the planet 
of Ag ~ 0.06 at small phase angles. Increasing the particle 
albedo to a = 0.99 dramatically increases the planet's geo- 
metric albedo to Ag = 0.55 at small phase angles. Particle 
albedos of a = 0.8 and a = 0.2 produce geometric albedos 
at zero phase oi Ag — 0.25 and 0.03 respectively. Our results 
are in agreement with the geometric albedos calculated by 
Dlugach & Yanovitskij ( 1974| who numerically solve for the 
intensity of radiation emerging from an atmosphere using 
iterative techniques. 

The single scattering phase function determines the an- 
gular distribution of scattered photons. In our initial inves- 
tigations we use the Henyey-Greenstein phase function de- 
scribed by the asymmetry parameter g. A strongly forward 
scattering phase function will scatter incident stellar pho- 
tons deeper into the atmosphere with the effect of reduc- 
ing the overall geometric albedo of the planet. Specifically 
we find zero phase geometric albedos can be reduced by 
a factor of three when the asymmetry parameter is varied 
from isotropic scattering {g = 0) to very forward scattering 
(5 = 0.9). 



3.2.2 Lower boundary condition 

The bottom layer of our simulation represents the bottom of 
the atmosphere. On a rocky planet, the bottom layer should 
act like a surface of a rocky planet, be it snow, water, or 
sand. On a gas giant planet, the bottom layer is only defined 



as the limit of our simulation. We must specify how much 
atmosphere to simulate, and how to treat the photon pack- 
ets that make it through the entire (simulated) atmosphere. 
We assign one of two characteristics to the bottom layer: 
Lambert or absorptive. A Lambert surface isotropically re- 
fiects all radiation, while an absorptive surface absorbs all 
the incident photons. Actual absorbed photons are often re- 
processed to longer wavelengths and do not contribute to 
the optical scattered light. 

Varying the bottom layer makes a huge impact on a 
small subset of our models with sparse cloud coverage or very 
low optical depths. In these situations the Lambert surface 
produces much larger geometric albedos than the absorb- 
ing lower boundary condition. For atmospheres with optical 
depths > 10 or porosity levels < 80%, the lower boundary 
does not affect our results since the vast majority of stel- 
lar photons interact with the upper levels of the atmosphere 
and do not penetrate to the bottom. 



3.3 Optical Depth effects in smooth and 3D 
atmospheres 

Figure [5] presents results of simulations for four different at- 
mospheric porosity levels and four different optical depths of 
the equivalent smooth model. Along with the smooth model, 
this figure therefore shows twenty different reflected light 
models covering what we believe to be a range of optical 
depth and porosity relevant to atmospheric cloud structures. 

The atmospheric optical depth determines photon pen- 
etration and is a major factor in determining a planet's ge- 
ometric albedo. With a small optical depth, photon packets 
pass through the atmosphere with no interaction and the 
lower boundary condition determines the reflectivity. 

As the optical depth increases, the 3D models and the 
smooth models react differently. The geometric albedos for 
smooth and 3D models initially increase with optical depth. 
Above T2 ~ 10 the smooth models maintain a constant geo- 
metric albedo. However, for 3D atmospheres the geometric 
albedo begins to decrease for high optical depths. This is a 
result of the increased depth of penetration of photons in a 
3D atmosphere compared to the equivalent smooth model. 
While photons that penetrate deep into the atmosphere may 
scatter back out via low optical depth paths, there is also 
a high probability that they will be absorbed by optically 
thick clumps encountered along their random path back out 
of the atmosphere. Recall from § |2.1| that the equivalent 
smooth model has the same total mass as a 3D model, but 
in the 3D model the mass is redistributed via the hierarchi- 
cal clumping algorithm resulting in regions of significantly 
lower and higher optical depths and thus an overall increase 
in depth of penetration. This is a well known result from 
three-dimensional radiation transfer (e.g., |Boisse| ( |1990| ). 

In addition to Tz , figure [5] shows that the atmospheric 
porosity has a significant impact on a planet's reflectivity. As 
described in § |2.7[ we vary the porosity of the 3D models and 
a high porosity results in only a few dense cloud structures in 
the grid. For lower porosities the hierarchically distributed 
mass results in a more continuous cloud structure and larger 
cloud covering factor (see Figure [2|. 

The main trends that we observe in Figure [5] are 1) 
greater cloud coverage results in a higher overall albedo, 2) 
compared to the equivalent smooth models, 3D atmospheres 
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Figure 5. Reflected light phase functions for our set of 3D atmospheric models with different porosity levels. All models here used 
a = 0.5, g = 0.5, absorbing bottom layer, /smooth = 0. Notice that 3D eff'ccts lead to a reduction in the reflectivity for all optical depths 
considered (the possible exception being r = 0.1, where the results of the smooth model are actually indistinguishable from several 
porous models, due to photon noise at that low level.) 



have a lower Ug, and 3) different porosities have different 
optical depths beyond which the geometric albedo ceases to 
rise and begins to decrease. 

Intuitively, if more of the atmosphere is covered with 
clouds, then the overall albedo will be higher. For an opti- 
cal depth r — 10, the simulation with 100% cloud coverage 
has a zero phase geometric albedo Ag — 0.08; while cloud 
coverages of 90%, 45%, and 15% give zero phase geometric 
albedos Ag = 0.06, 0.03, and 0.01 respectively. This is ex- 



pected because more cloud coverage will result in more pho- 
tons being reflected from the upper layers of the atmosphere. 
Note these simulations have an opaque lower boundary and 
/smooth = 0. As /smooth IS iucreascd the geometric albedo 
increases towards the values for the smooth models. 

When compared to the equivalent smooth models, all 
our 3D simulations with significant optical depth and inter- 
nal structure have lower geometric albedos. An albedo bias 
between uniform and one-dimensionally fractal models has 
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been reported in atmospheric science literature, specifically 



modeling the Earth's atmosphere (Cahalan 1994 Cahalan 
et al.|1994 |. If treated as plane-parallel, California stratocu- 
mulus clouds have a relative albedo bias of 15%, meaning 
that smooth simulations will suggest an albedo 15% higher 
than the albedo of the actual clouds, which have a great deal 
of fractal structure. 



3.4 Discussion of Scattered Light Levels 

Our investigation finds that 3D radiation transfer effects in 
an inhomogeneous atmosphere lead to lower levels of re- 
flected light compared to uniform atmosphere models of the 
same total mass. For inhomogeneous atmospheres low den- 
sity paths allow the photons to penetrate deep into the at- 
mosphere where there is a higher probability they will be 
absorbed in a dense cloud on their path back up through 
the atmosphere. This lowers the geometric albedo and is 
a 3D radiation transfer effect. However, perhaps the most 
important single parameter is the albedo of the individual 
scattering particles in the atmosphere (Sudarsky et al.|2005 



Barman et al. 2005 Seager et al. 2005 1. Very high albedos as- 



sociated with enstatite and iron particles (typically a > 0.9, 
see next sections) can result in reflected light levels an order 
of magnitude higher than an atmosphere where the particles 
have, for example, a — 0.5. 

Overall, the geometric albedos are much lower than 
expected and make reflected light detection very difficult. 
Other papers have specifically discussed the low geomet- 
ric albedos of extrasolar planetary atmospheres, including 



Marley et al. (19991, Seager et al. (20001, while Dlugach & 



Yanovitskij|(|1974[) found low theoretical albedos for generic 



planets. Marley mentions that, using their model, planets 
warmer than 400 K will be dark and that most planets will 
be dark in reflected light beyond a 600 nm wavelength. Sea- 
ger shows that at wavelength 480 nm the geometric albedo 
of her model varies between below 0.01 and 0.4, depending 
on the mean particle size. 



Finally, Dlugach & Yanovitskij ( 1974 1 used an analyt- 



ical method of calculating the intensity of radiation from a 
homogeneous atmosphere to flnd that isotropic scattering in 
an atmosphere with particles with a single scattering albedo 
of 0.9 would only result in a geometric albedo of 0.3. Their 
models speciflcally illustrate the necessity of a very high sin- 
gle scattering albedo to produce a high geometric albedo. In 
a forward throwing atmosphere, a seemingly minute change 
in the single scattering albedo from 0.99 to 0.98 can change 



the geometric albedo from 0.34 to 0.26 (Dlugach & Yanovit- 
skij|1974 |. It is possible that many researchers were expect- 
ing (or hoping) extrasolar giant planets would have the same 
high geometric albedos as Jupiter {Ag — 0.52) and Saturn 
{Ag = 0.47). However, at the temperatures of the extraso- 
lar planets, there seem to be precious few condensates with 
the high single scattering albedos. This is a sharp contrast 
to the condensates of the gas giants of our own solar sys- 
tem, which typically have single scattering albedos of 0.996 
( Chanover et al.||1997 |. Much as we were expecting to flnd 
Jupiters around other stars, and instead found close-in gi- 
ant planets, we were expecting to flnd Jupiter-like geometric 
albedos and instead are flnding dark, low geometric albedos. 



4 REFLECTED LIGHT FROM HD 209458B 

In the previous section we explained how inhomogeneous 
planetary atmospheres lead to very low geometric albedos 
compared to smooth atmospheres. In this section we con- 
struct a 3D scattered light model for HD 209458b, arguably 



the best studied extrasolar planetary atmosphere (Seager 
fc Sasselovl[2000| |Seager|[20M| [Deming et aL]|2005| |Sho^ 



man et al. 



20081. The goal of this model is to investigate 



3D radiation transfer effects on the reflected light levels, in 
particular the geometric albedo limit from the MOST satel- 
lite. In addition, using recent ID atmospheric models, we 
investigate limits on the atmospheric parameters such as 
cloud height, location, and composition. As with the mod- 
els presented earlier, we do not construct a self-consistent 
3D temperature-pressure-density model for HD209458b, but 
take a ID atmosphere model and convert it to 3D using 
the hierarchical clumping algorithm described earlier. There 
have been recent 3D atmosphere models for HD209458b, but 
these models did not include 3D radiation transport. There- 
fore, what we present in this paper are the flrst 3D scattered 
light models for extrasolar planets and may be seen as a flrst 
step towards incorporating accurate 3D radiation transfer 
into global 3D atmospheric models. 

For our HD209458b models we extended the Monte 
Carlo simulations to include exponential density gradients, 
multiple opacity sources mixed throughout the simulation 
grid, and scattering from pre-tabulated phase functions as 
may be computed from Mie scattering codes. In addition, 
we can limit the lower and upper boundaries of a cloud of 
scattering particles to specific heights in the atmosphere. 

In this section we detail how we determined our fiducial 
choices for the parameters for individual scattering species. 
We explain which parameters might reasonably be varied, 
and we simulate those variations to compare the geomet- 
ric albedo of the varied atmospheres. We present results for 
scattered light at 500nm, at which wavelength absorbed pho- 
tons will not primarily contribute to the optical fiux. 



4.1 HD209458b atmospheric parameters 

We must first compute the specific parameters of the simula- 
tions, such as the densities and opacities of the condensates. 
Two of the most important condensates in HD 209458b's 
atmosphere are believed to be enstatite (MgSiOa), with an 
albedo of 0.999 and iron (Fe), with an albedo of 0.685 [AcF| 



erman & Marley (20011. They both condense high in the 



atmosphere, which makes them significant contributors to 
the geometric albedo. Further, the lack of other high albedo 
species (in the temperature's and pressure's of HD 209458b) 
indicates that enstatite must be present in very reflective 
planets. Our focus in simulating HD 209458b will be on sim- 
ulating a Rayleigh scattering atmosphere mass-dominated 
by hydrogen, with Mie-scattering clouds of both enstatite 
and iron. In all cases we will be simulating photon packets 
at 500 nm, and choose the optical properties (opacity, scat- 
tering, and phase function) accordingly. This is a reasonable 
choice of simulated wavelength due to HD 209458 being a 
GOV star. 
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4.1.1 Density profile 

The smooth density structure we initiaUy adopt is pro- 
vided by the ID atmosphere code described by |Seager et al.| 
( 2005 1 which self consistently solves the equations of radia- 



tion transfer, radiative equilibrium, and hydrostatic equilib- 
rium. (We should mention that there is some indication that 
HD 209458b is not in hydrostatic equilibrium at the top of 
the atmosphere, and this is a limitation of our simulations 
I Vidal-Madjar et al.|2004 |.) Using the temperature-pressure 
profile from Seager's code we compute the vertical hydrogen 
density structure assuming the ideal gas law holds through- 
out the atmosphere. It is this ID density profile that we will 
convert to 3D using our hierarchical clumping algorithm. 

We find the mass density profile of hydrogen computed 
from Seager's code fits an exponential of the form. 



3.4 X lO'^^e'^/^gcm"^ 



(3) 



where p is the density, z is the height in cm in the atmo- 
sphere, and the scaleheight H = 1.04 x lO^cm. The total 
height of the atmospheric simulation is determined from hy- 
drostatic equilibrium to be 8.2 x lO^cm, about 9% of the 
radius of the planet. 

Assuming solar abundances for the elements Mg, Si, 
and O, we construct the smooth ID density structures of 
enstatite and iron. For every 10^^ atoms of hydrogen, there 
are 10*-^^ atoms of O, 10^-^* atoms of Mg, and lO'''^'^ atoms 



of Si, so the Si is our limiting factor ( Cox 2000 1. The density, 
p, which we computed above is the density of H2 , which has 
a molecular weight of 2.02 gmole"^. 

Once we have the moles of H, we can determine the 
moles of MgSiOs, since using the above ratios, we know that 
for every one H atom we can have 3.548 x 10~^ molecules of 
MgSiOa. We also have a free parameter, c, which we will use 
to describe the percentage of enstatite that has condensed 
into the cloud. Thus, our density structure of enstatite con- 
densed in the atmosphere is: 



PMgSiOa = c X pH X 3.59 X 10 



(4) 



Similarly, we calculate the density of iron in the atmo- 
sphere assuming solar abundanc es. For ever y lO'^^ atoms of 
H, there are lO''-^" atoms of Fe ( |Cox||2000| ). Using this ra- 
tio, for every H atom we will have 3.467 x 10~^ atoms of Fe 
giving. 



pFc — cx ph X 1.917 X 10 



(5) 



where again c is the condensation fraction. 

As noted in § |2.2[ we previously made individual simu- 
lations in cubes, then we tiled those cubes onto the surface 
of a sphere in order to represent the entire planet. However, 
in order to match the height of the simulated data, we have 
had to remake our individual simulations into a rectangular 
grid, with 2 being five times the length of x or y. We main- 
tained the cubic dimensions of our individual grid cells, and 
instead simply use five times as many grid cells in the z 
direction (500 z cells instead of 100 cells in x and y). The 
height of our entire atmospheric simulation is 8.235 x 10* 
cm, and our individual grid cells are 1.647 x 10® cm on each 
side. Each of the cubic cells has a volume of 4.467 x lO^^cm^. 



4.1.2 Location of enstatite and iron clouds in atmosphere 

To find the placement of the enstatite cloud in our simulated 
atmosphere, we used the condensation curves from figure 4 
of |Seager et al.| ( |2000[ ), which suggest that in our tempera- 
ture pressure profile, enstatite condenses at approximately 
4.1 X 10* cm, or about halfway up the atmosphere. We chose 
to have the enstatite cloud continue through two pressure 
scale heights, similar to Jupiter's clouds, which places the 
upper boundary of the cloud at approximately 5.8 x 10* cm, 
or 70% the height of the atmosphere. 

Using the same techniques as for enstatite, we find that 
the iron should condense at 3.4 x 10* cm, or about 42% 
up the atmosphere. We again chose to have the iron cloud 
continue through two pressure scale heights, which puts the 
upper boundary of the iron cloud at approximately 5.1 x 10* 
cm, or 62% the height of the atmosphere. 

The clouds we simulate are fractal, largely without 
sharp edges, and they follow the same pressure profile, mean- 
ing they will usually be optically thicker at the bottom than 
at the top. 



4.1.3 Opacity and scattering properties of hydrogen, 
enstatite, and iron 

Hydrogen opacity is calculated from the cross section calcu- 
lated by Dalgarno fc WiUiarns| ( |1962| |, 



a(A) = 



8.14 X 10" 



1.28 X 10" 
AS 



+ 



1.61 
AS 



(6) 



where the wavelength A is in A and the cross section is in 
cm'^. Further, using the relationship an — up where n is 
the number density, we find the opacity of molecular hy- 
drogen at 500 nm is k = 4.14 x lO^^cm^g"^. Knowing the 
density and opacity of molecular hydrogen, and assuming 
Rayleigh scattering, allows us to include molecular hydrogen 
in our simulations. Though we generally assume hydrogen 
acts only as a scatterer (an albedo of 1), we assigned hydro- 
gen an albedo of 0.999, so that photons that passed through 
the enstatite and iron clouds would not scatter around the 
optically thick hydrogen bottom layer indefinitely, but in- 
stead would eventually be absorbed. This is a reasonable 
approximation because it is unlikely that any atmosphere 
will consist of pure hydrogen and will likely contain absorb- 
ing contaminants. 

The opacity and albedo of enstatite and iron are cal- 



culated using a Mie scattering code (Bohren & Huffman 



19831. For enstatite at 500nm the optical constants n and 
k are 1.577 and 2.1 x 10~^, respectively (Dorschner et al. 



1995 1 . Using these constants and assuming the enstatite 



is in the form of 5 pm diameter spheres, we calculate the 
opacity of enstatite to be 1.12 x 10*cm^g~^ and the albedo 
"MgSiOa = 0.999. Using the optical constants of iron (jOrdal 
et al.|1985l, n = 2.74 and k = 2.88, along with the assump- 



tion that the condensates are 5 pm in diameter, allows us to 
use Mie theory to calculate the optical parameters of our sin- 



gle sc attering iron particles at 500 nm ( Bohren & Huffman 
1983 j ). We find that the opacity of iron is 4.03 x lO^cm^g"^ 
and that the albedo is apc = 0.685. 

The scattering phase functions for enstatite and iron 
are calculated from Mie theory and tabulated for the Monte 
Carlo scattering code. Using pre-computed tabulated phase 
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Figure 6. Henyey Greenstein versus Mie scattering phase func- 
tions for iron and enstatite. The Mie phase functions are calcu- 
lated for five micron radii enstatite and iron grains. The enstatite 
has a slight back scatter around 180 degrees. The Henyey Green- 
stein g=0.9 phase is similar to the enstatite phase, but has no back 
scatter. The g=0.9 phase also underestimates the iron phase at 
180. We opted to use the Mie phase functions in this simulation 
of HD 209458b. 



functions allows a wide variety of phase functions to be in- 
corporated into our scattering codes, and for enstatite allows 
us to have a slight back-scattering component not possible 
with a Henyey-Greenstein phase function (see figure [m. 



4.2 Results 

Before reporting our results, we should note that our fidu- 
cial ID model is the only self-consistent atmospheric struc- 
ture model we use. In our 3D scattered light simulations, we 
freely move the condensate clouds, adjust the opacities (by 
changing the mass), and change the absorption of the gas. 
We do this primarily because we want to isolate the effects 
of our changes. We also do this because there does not yet 
exist a fully self-consistent 3D atmosphere model that accu- 
rately includes 3D radiation transfer effects, although there 
are moves to address this( Showman et al.|2008 l. Our varia- 
tions can be equated to small changes in the modeled atmo- 
sphere. For instance, changing the height of the clouds can 
be equated to using a different temperature pressure profile. 
Changing the absorption of the gas assumes that a species 
other than molecular hydrogen might be prevalent in the at- 
mosphere. Using a different albedo for a condensate tests the 
effects of an impure (or altogether different) condensate in 
the atmosphere. By changing these individual input param- 
eters, we no longer create a self-consistent model, but we are 
able to isolate the effect of our variations. We feel that with 
this first 3D modeling of scattered light from an extrasolar 
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Figure 7. Density X opacity schematic of the atmospheric sim- 
ulation for HD 209458b. This is a column summation for all y. 
Our entire simulation is approximately 8000 km in height, with 
the temperature and pressures noted at the left. The dark re- 
gions of the schematic are optically thickest while the white re- 
gions are optically thin. Molecular hydrogen is throughout the 
exponential density atmosphere. Enstatite clouds condense high- 
est in the atmosphere, and overlap with the iron which condenses 
slightly below them. Certain features have been exaggerated in 
this schematic, such as the opacity due to hydrogen, and the tenu- 
ity of the enstatite and iron clouds, in order to see the features 
more clearly. 



planetary atmosphere, isolating the effect of our variations 
is more important than a self-consistent simulation. 



4.2.1 Fiducial 3D model 

Using the parameters computed above for the ID den- 
sity structure, we have created a fiducial 3D model of HD 
209458b which generates a geometric albedo Ag = 0.42 at 
500 nm. The atmospheric scattered light simulation con- 
sists of three scatterers: Rayleigh scattering hydrogen gas, 
Mie scattering five /xm enstatite, and Mie scattering five /xm 
iron. The pressure and temperature of HD209458b are pro- 
vided by Seager's ID calculation from which the density is 
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Figure 8. Geometric albedo phase curve of our fiducial model, 
including a two pressure scale height enstatite cloud and a similar, 
but lower, iron cloud in the atmosphere, with Rayleigh scattering 
hydrogen gas throughout. 



Figure 9. Geometric albedo phase curve of varying abundance 
rates for the iron and enstatite clouds. We vary some values above 
100% because HD 209458b may be more metal-rich than our sun, 
resulting in more enstatite in the cloud. 



computed. We assume a 10% condensation rate, c, for both 
enstatite and iron. We convert this ID model to 3D as de- 
scribed above, assume a porosity (as described in § 2.2 I of 
45%, and we assume that the particulates will condense ac- 
cording to their condensation curves and will extend up two 
pressure scale heights (Seager et al. 20051. Our 3D radia- 



tion transfer simulation models the outermost ~ 8000 km of 
the atmosphere. It is against this fiducial model (see figure 
[7| that we will compare our other simulations. With these 
fiducial parameters, the resultant geometric albedo of HD 
209458b is Ag ~ 0.42, as seen in figure [8] 



4-2.2 Abundances and condensation 

To determine the extent of condensed enstatite and iron 
in the atmosphere, we use a single abundance parameter 
that encompasses both condensation and abundance. For 
instance, in our model, an atmosphere with lOx solar abun- 
dance and 10% condensation is treated the same as a Ix 
solar abundance with a 100% condensation. 

For our fiducial model, we assumed that 10% of the 
enstatite and iron in the atmosphere would condense into 5 
/im grains. The condensation of enstatite and iron could be 
affected by the purity of the metals in the atmosphere, as 
well as the turbulence. 

Our fiducial model initially assumes a solar abundance 
of the constituents of enstatite and iron. However, the plan- 
etary system as a whole could have a higher abundance 
than solar, as HD 209458 has a higher abundance than so- 
lar ( |Schneider|[2006[ ). Additionally, HD 209458b could have 
a higher or lower abundance of enstatite or iron than solar, 
as Jupiter and Saturn have some metals at differing concen- 



condensates could form at higher regions in the atmosphere 
and rain down into the cloud deck, creating locally higher 
or lower concentrations of enstatite and iron than the planet 



as a whole Ackerman & Marley ( 2001 1 



tration rates Guillot (19991. Finally, it is also possible that 



Because our model only considers this single parame- 
ter which encompasses both condensation and abundance, 
we vary our enstatite and iron levels widely to encompass 
both low abundances and condensation rates, as well as high 
abundances and high condensation rates. We vary this pa- 
rameter from 0.001 to 10 times solar abundance. 

Figure |9] demonstrates that the abundance of enstatite 
and iron in the cloud only has a small impact on the geo- 
metric albedo, when it is above a minimum amount, around 
1% of solar abundances. Similar to our earlier, more gen- 
eral, models (see § Isl, larger opacities combined with a 
fractal cloud lower the overall reflectivity. This leads to a 
100% abundance producing a slightly less reflective geo- 
metric albedo than our 10% fiducial model. This is due to 
the light entering the cloud at optically thin regions but 
not escaping because it is surrounded by optically thick re- 
gions. However, even the 1000% abundance model shows 
only modest differences in the geometric albedo from the 
fiducial model, within the noise limits of randomly changing 
the cloud structure but maintaining the same porosity. 



4.2.3 Porosity 

Our fiducial scattered light model assumed that the cloud 
porosity is 45% and we have compared this with models with 
porosities of 80% and 2%. It should be noted that when 
we consider percent porosity, we are only considering the 
grid cells containing the condensate cloud, the hydrogen is 
considered to be smoothly distributed with the exponential 
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Figure 10. Geometric albedo phase curve of varying porosity 
rates. The nearly smooth model (2% porosity) returns a geomet- 
ric albedo only slightly better than our fiducial model of 45% 
porosity. The more tenuous 80% porosity model gives a much less 
reflective planet. 



density described above. In all of our porous models the 



Figure 11. Geometric albedo phase curve of cloud heights. Our 
fiducial (predicted) model puts the enstatite cloud highest in the 
atmosphere, with the iron below it. If we instead mix the two 
clouds, and put both iron and enstatite at the same height in 
the atmosphere, the geometric albedo is reduced. In the situation 
where we switch the placement of iron and enstatite, with iron 
on top, the geometric albedo is significantly reduced, to only 10% 
the value of our fiducial model. 



As expected from our general investigations, figure [To] 
shows that cloud porosity is a significant factor in determin- 
ing the overall reflectivity of HD209458b. If incident starlight 
is able to travel deeply into a highly porous cloud, through 
the optically thin regions, then that light is often not able 
to emerge from the atmosphere, and results in a lower geo- 
metric albedo. The porosity is even more important in the 
specific situation of HD 209458b than it is in our general 
cases in § [3] This is because of the relatively low albedo of 
iron (ape = 0.685) compared to enstatite (ajvigSiOa = 0.999). 
In general the enstatite cloud will condense at a higher al- 
titude than the iron cloud, meaning that in a smooth cloud 
deck the photons will mainly interact with the optically thick 
enstatite and never reach the lower iron cloud. However, 
a porous enstatite cloud will allow starlight to penetrate 
deeper into the atmosphere, where it may be absorbed by 
the low albedo iron. By comparison, in our general investi- 
gations changing the porosity from 45% to 80% reduced the 
geometric albedo by about 30%. The same change in the HD 
209458b atmosphere reduced the geometric albedo by 40%. 



4.2.4 Cloud height 

The height at which the enstatite and iron clouds form may 
be varied and this is equivalent to exploring situations in 
which the atmospheric clouds have not settled as well as 
predicted from our ID temperature-pressure model. Cloud 
heights may differ due to a very turbulent atmosphere, or 
because there is a higher abundance of iron than we suppose. 
Additionally, changing the cloud height is a way we can test 



different temperature pressure profiles. We are aware that 
this is not self-consistent, in that we are placing the clouds 
at temperatures and pressures where they would not oc- 
cur, but this still gives us good first approximations as to 
what would happen with completely different parameters. 
We tested cases in which the iron and the enstatite occur at 
the same cloud height, and in the unlikely situation of iron 
occurring at the enstatite height in the atmosphere, and en- 
statite occurring at the iron height. 

Figure 111] shows the height of the condensate clouds 
hugely affects the geometric albedo of the planet, for much 
the same reason as the porosity rate affects it. In this case, 
it is mainly the height of the iron cloud that makes the 
difference. Its low albedo can sharply decrease the overall 
reffectivity of a planet, producing an albedo only 14% as 
reffective as the fiducial model. 



4.2.5 Gas absorption 

To simulate additional absorptive constituents in the hydro- 
gen gas component of the model atmosphere we ran scat- 
tered light simulations in which the gas albedo is lowered 
from 0.999 to 0.5. This is a possibility if TiO, CH4, or Na 
molecules are in the atmosphere, because they are much 
more absorptive than hydrogen. The simulations also ex- 
plored the effects on the planet's geometric albedo of raising 
and lowering the height of the enstatite cloud. 

Figure [12] shows there is very little difference between 
the case where the enstatite clouds were high in the atmo- 
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Figure 12. Geometric albedo phase curve of different absorption 
for the gas, and different heights of the enstatite cloud. Chang- 
ing only the gas absorption does not affect the geometric albedo 
much, nor does making the enstatite cloud lower. However, if we 
make the gas more absorptive and lower the cloud, then the light 
is absorbed by the gas before it reaches the enstatite cloud, lead- 
ing to a lower geometric albedo. 



Figure 13. Geometric albedo phase curve of different single scat- 
tering albedos for enstatite. The relatively small change in the 
single scattering albedo, from 0.999 to 0.985 has a large impact 
on the geometric albedo of the planet. This change in single scat- 
tering albedo could result from more amorphous than glassy en- 
statite in the atmosphere. Note that we kept the same phase 
function for enstatite, and only changed the albedo. 



sphere and the gas albedo was varied. The opacity of the 
molecular hydrogen is too small to make much of a differ- 
ence in the geometric albedo when most photons instead in- 
teract with the much more opaque and higher enstatite and 
iron. However, the case changes when the clouds are moved 
to the bottom of the atmosphere. Even though the clouds 
are much more opaque than the hydrogen, the photons still 
interact and are absorbed if the cloud is around 2% of the 
height of the atmosphere instead of 52% high. Merely low- 
ering the cloud or lowering the molecular hydrogen albedo 
will not change the geometric albedo, but the two changes 
together cut the geometric albedo by more than half. 



times, a small change in the albedo contributes significantly 
to the geometric albedo, as we show in figure [131 Comparing 
glassy and amorphous enstatite is outside the scope of this 
work, but it is sufficient to say that the constituents of ex- 
trasolar planetary atmospheres are not known well-enough 
to favour one or the other. It should be noted that though 
we used the adjusted optical constants to compute the sec- 
ond albedo, we did not recompute the phase function of the 
enstatite based on the new constants. This was to isolate the 
change in albedo and we see that what appears to be a minor 
albedo difference, between amorphous and glassy enstatite, 
can nearly halve the geometric albedo of the planet. 



4-2.6 Enstatite albedo 

By using a slightly lower value for the single scattering 
albedo of enstatite, our models produce a geometric albedo 
of only Ag = 0.25. Using Mie theory and optical constants 



for glassy enstatite from Dorschner et al. ( 1995 1, we calculate 
a single scattering albedo for enstatite of aMgSiOa = 0.999. 



However, optical constants by Scott fc Duley| ( 1996 1 dif- 



fer slightly, and measure amorphous enstatite. These val- 
ues are used by [Marley et al.| ( |1999| to develop their ex- 
trasolar giant planet atmospheric model. Using the values 



of Scott & Duley ( 1996 1 results in an albedo of enstatite 



of aMgSiOa = 0.985. This minute change contributes signifi- 
cantly to the geometric albedo, because the phase function 
of enstatite throws most photons forward, so most photons 
penetrate deep into the atmosphere and scatter many times 
before coming out. Because the photons scatter so many 



4.3 Discussion of HD 209458b Model Atmosphere 

For HD 2Q9458b to have a high geometric albedo, the 
enstatite cloud condensate must be high in the atmo- 
sphere, smooth, and pure. At the temperatures and pres- 
sures present in the atmosphere of HD 209458b, there are 
very few condensates which can make a positive impact on 
the geometric albedo of the planet ( Seager et al.|2000 l. Ac- 
cording to the condensation curves, enstatite forms high 
in the atmosphere and has a single-scattering albedo high 
enough to give the planet a high geometric albedo. 

However, the geometric albedo of HD 209458b has a Sc- 
upper limit of 0.17 in the visible wavelengths, as determined 
by the MOST sateUite ( |Rowe et al.|2006[ 12007] ). This upper 
limit on the geometric albedo seriously constrains the na- 
ture of the enstatite in the atmosphere of HD 209458b. The 
enstatite cannot be high, smooth and pure as we speculated 
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in our fiducial model, because if it were, then the reflected 
starlight would have been detected by MOST. 

We find several variations of our models which depress 
the geometric albedo of HD 209458b. Lowering the abun- 
dance (comprised of both condensation and abundance com- 
pared to solar) from 10% to 0.1% cuts the geometric albedo 
by more than one half. Very porous clouds lower the albedo 
by one quarter. Significantly lowering the enstatite cloud 
and raising the gas absorption together more than halve the 
geometric albedo, though each change on its own does not 
affect the geometric albedo. Lowering the single scattering 
albedo of enstatite to 0.985 from 0.999 reduces the geomet- 
ric albedo Ag = 0.25 compared to our fiducial model value 
Ag = 0.42. Finally, mixing the high-albedo enstatite with 
the low-albedo iron more than halves the geometric albedo, 
and if iron is higher in the atmosphere, then the geometric 
albedo is only Ag — 0.04. 

Though a 0.1% abundance of enstatite results in a low 
geometric albedo for the planet, we are unconvinced that a 
low abundance of enstatite is what is causing the low albedo 
of HD 209458b. First, the metallicity of HD 209458 is higher 
than solar metallicity, suggesting that the enstatite in the 
system could be higher than solar ( Schneider|2006 l. Second, 
we notice the enhanced metallicity in the atmospheres of our 
own gas giant planets, Jupiter and Saturn, compared to so- 
lar, and consider it reasonable that extrasolar giant planets 
will also exhibit an enhanced metallicity compared to the 
parent star. Finally, we think the cloud condensation can 
promote a locally enhanced abundance at the condensation 
height of enstatite in the atmosphere. 

If clouds are very porous, the geometric albedo is re- 
duced, but not enough to explain the low geometric albedo. 
In our three test simulations, the most porous case with 80% 
porosity still gives a geometric albedo Ag = 0.3, above the 
upper limit of the MOST data. 

If the cloud is much lower in the atmosphere than the 
condensation curves suggest, and the atmosphere is filled 
with more absorptive Rayleigh scatterers, then the geomet- 
ric albedo is reduced to Ag ~ 0.18, still above the upper 
limit imposed by MOST. This is a reasonable situation if 
ions or other species expected in the atmosphere are highly 
absorptive. However, as mentioned before, the opacities of 
the condensates are much higher than the opacities of the 
Rayliegh scatterers, so the condensates must be very low in 
the atmosphere, and the Rayleigh scatterers very absorptive 
to generate a low geometric albedo. 

The nature of the enstatite, either glassy or amor- 
phous, significantly affects the geometric albedos, because 
it changes the single scattering albedo from 0.999 to 0.985. 
This small change reduces the geometric albedo from Ag — 
0.42 to 0.26, above the geometric albedo limit of MOST. We 
do find it hard to believe that all the enstatite condensed in 
the atmosphere is pure, glassy enstatite. Other species are 
present in the atmosphere, even some other condensates, so 
it is a reasonable expectation that some of the reduction 
in geometric albedo comes from amorphous and impure en- 
statite depressing the single scattering albedo of the conden- 
sate, but obviously, just the change to amorphous enstatite 
is not sufficient, because the geometric albedo remains above 
our upper limit. 

The atmospheres of Jupiter and Saturn are convective 
and as a result mostly homogenous. The upper layers are 



not completely homogenous, due partly to condensation of 
clouds, but there is still diffusion of heavy elements through- 
out the atmosphere ( Guillot|1999 1 . In our simulations where 
the iron and the enstatite have mixed, the geometric albedo 
drops to less than half the value of our fiducial model. If 
iron is above the enstatite, the geometric albedo drops to 
Ag — 0.04. Several groups comment on the turbulence ex- 



pected in the atmosphere of HD 209458b ( Showman & Guil- 
lot|2002[ [Cho et aH2003J . We think that a turbulent atmo- 
sphere, with routine mixing of various condensates (and the 
resulting drop in the average single scattering albedo), is a 
significant contributor to the low geometric albedo. 



5 SUMMARY 

The primary goal of this paper was to investigate the effects 
of 3D radiation transfer on the optical reflected light levels 
from extrasolar planetary atmospheres. We find that com- 
pared to smooth models, 3D atmospheres reflect less light 
because starlight can penetrate deep into the atmosphere, 
increasing its probability of being absorbed. As found in 
other ID and 2D scattered light simulations, other atmo- 
spheric parameters serve to reduce the geometric albedos of 
extrasolar planets, most notably the scattering albedo of in- 



dividual particles in the planetary atmosphere (Dlugach & 



Yanovitskij||1974 \ . Indeed, we suspect that it is the albedo 



of particles in the atmosphere that results in the very low 
geometric albedo observed from HD209458b. 

We realize that our model is not a perfect representa- 
tion of an extrasolar atmosphere, and that other effects not 
explored here might contribute to a lower geometric albedo. 
Specifically, strong absorption of alkali metals sodium and 
potassium could lower the albedo throughout the visible 
spectrum (Sudarsky et al. (20001 and Soagcr et al. (2000)). 



However, we have tried to isolate some effects, specifically 
from the 3d structure of the atmosphere, that might con- 
tribute to the overall albedo of extrasolar planets. 

HD 209458b is a dark planet, with a geometric albedo 
Ag < 0.17. One possibility is that clouds are absent and 
the low albedo is due to atomic or molecular absorption. To 
study the albedos in the presence of clouds, we investigated 
a fiducial 3D scattered light model which considered the 
densities, opacities, and albedos of hydrogen, iron, and en- 
statite, and their respective placements in the atmosphere. 
This fiducial model gives a geometric albedo Ag > 0.4, so we 
know that it does not accurately represent the atmosphere 
of HD 209458b. Upon varying parameters in our model we 
found significant drops in the geometric albedo. Parameters 
that strongly affected the geometric albedo include mak- 
ing the clouds more porous, lowering the condensation rate 
of enstatite, and assuming amorphous enstatite instead of 
glassy. Indeed, we were surprised by how easily the geomet- 
ric albedo could be reduced with a relatively small change 
to the constituents of the atmosphere. However, we feel the 
most important parameter influencing the reflected light lev- 
els is the purity of enstatite and the assumption that it sits 
alone and above all the other condensates in the atmosphere. 
Models have suggested that HD 209458b has a turbulent at- 
mosphere, and if clouds are indeed present we expect that 
mixing of enstatite with other, lower albedo condensates is 
what keeps the geometric albedo so low. 
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